cycling/ros_variance.m

   1 function [ros,ros_mean] = ros_variance(n,ps)
   2 %estimates variance in ros between two detections
   3 %n - number of random points to use in simulation
   4 %ps - struct with path information, ps = cluster_paths(w,1)
   5
   6 real_data = input_num('Use real data? 0 = no 1 = yes',1,1);
   7
   8 if real_data == 0
   9     %coordinates of detections
  10     d1 = [0 0 0];
  11     d2 = [3 4 12];
  12
  13     %reported
  14     dist_r = sqrt((d2(1)-d1(1))^2+(d2(2)-d1(2))^2);
  15     time_r = d2(3) - d1(3);
  16     ros = dist_r/time_r;
  17     fprintf('ROS = %f\n',ros)
  18
  19     %statistics
  20     sig1 = 1;
  21     sig2 = 1;
  22     time_diff = 0;
  23
  24 else
  25     %real data
  26     %load p_struct.mat
  27     %ps = p;
  28     %clear p
  29     num_paths = length(ps.paths);
  30     path_num = max(2,round(num_paths*rand));
  31     p = ps.paths(path_num);
  32     path_length = length(p);
  33     path_end = max(2,round(path_length*rand));
  34     p1 = p.p(path_end -1);
  35     p2 = p.p(path_end);
  36     fprintf('Generating stats for path %d, points %d and %d \n',path_num,path_end-1,path_end);
  37     lon1 = ps.points(p1,2);
  38     lat1 = ps.points(p1,1);
  39     time1 = ps.points(p1,3);
  40     lon2 = ps.points(p2,2);
  41     lat2 = ps.points(p2,1);
  42     time2 = ps.points(p2,3);
  43
  44     E = wgs84Ellipsoid;
  45     dist_r = distance(lat1,lon1,lat2,lon2,E);
  46     time_r = time2-time1;
  47
  48     %testing numbers
  49     dist_r = 1500;
  50     time_r = 0.9;
  51
  52     ros = dist_r/(time_r*24*3600);
  53     fprintf('Dist = %f m Time = %f s \n',dist_r,time_r*24*3600);
  54     fprintf('ROS = %f m/s \n',ros)
  55
  56     %coordinates of detections
  57     d1 = [0 0 0];
  58     d2 = [0 dist_r time_r];
  59
  60     %statistics
  61     sig1 = 375;
  62     sig2 = 375;
  63     time_diff = 0.25;
  64 end %if real_data == 0
  65
  66
  67
  68 %n = 200;
  69 r1 = sig1*randn(n,3)+d1;
  70 r1(:,3) = d1(3)-time_diff*rand(n,1);
  71 r2 = sig2*randn(n,3)+d2;
  72 r2(:,3) = d2(3)-time_diff*rand(n,1);
  73
  74 dr = r2-r1;
  75 dist_rand = sqrt(dr(:,1).^2+dr(:,2).^2);
  76 time_rand = dr(:,3);
  77 if real_data == 0
  78     ros_rand = dist_rand./time_rand;
  79 else
  80     ros_rand = dist_rand./(24*3600*time_rand);
  81 end
  82 if n < 1000
  83     figure,scatter3(r1(:,1),r1(:,2),r1(:,3),'*b')
  84     xlabel('x'),ylabel('y'),zlabel('time')
  85     title('Random paths between two detections')
  86     hold on
  87     scatter3(r2(:,1),r2(:,2),r2(:,3),'*r')
  88     for i = 1:n
  89         plot3([r1(i,1) r2(i,1)],[r1(i,2) r2(i,2)],[r1(i,3) r2(i,3)],':')
  90     end
  91     hold off
  92     %make joint pdf
  93     figure,scatter3(dist_rand,1./time_rand,ros_rand);
  94     xlabel('Distance'),ylabel('1/Time'),zlabel('ROS')
  95 end
  96
  97 figure,histogram(ros_rand)
  98 if real_data == 0
  99     title('Histogram of ROS')
 100 else
 101     t_str = sprintf('Histogram of ROS \n Path: %d Points: %d and %d',path_num,path_end-1,path_end);
 102     t_str2 = sprintf('Distance: %f [m] Time: %f [h]',dist_r,time_r*24)
 103     title({t_str,t_str2})
 104 end
 105 xlabel('ROS'),ylabel('Number')
 106 ros_mean = mean(ros_rand);
 107 fprintf('Mean random ROS: %f m/s \n',ros_mean)
 108 ros_std = std(ros_rand);
 109 fprintf('Standard deviation of random ROS: %f m/s \n',ros_std)
 110
 111 figure,qqplot(ros_rand)
 112 title('QQ-plot')
 113
 114 end