detection/plot_like_new.m

   1 %function plot_like_new
   2
   3         % function h=heat(T)
   4         c = 2;  % time constant (h) : how long the heat takes to drop to exp(-1)
   5         %********************************
   6         heat = @(T) (T<=0).*exp(min(T,0)/c);
   7         %********************************
   8         % T = fire arrival time
   9         % 0 if T > T_now = 0, 1 at T=0, exp decay for T>= 0
  10         figure(1);fplot(@(T)heat(-T),[-c,5*c]);
  11         xlabel('time since fire arrival (h)');ylabel('heat fraction (1)')
  12         title('Burn model')
  13         fig
  14
  15         % function p=prob(h)
  16         false_detection_rate = 0.001;  % = 1/(1+exp(b))
  17         b = log(1/false_detection_rate  -1);
  18         half_prob_heat = 0.01; % 1+exp(-a*h + b) = 2 for h=half_prob_heat
  19         a = b / half_prob_heat;
  20         %********************************
  21         prob = @(h) 1./(1+exp(-a*h + b));
  22         %********************************
  23         err1 = false_detection_rate - prob(0)
  24         err2 = 0.5 - prob(half_prob_heat)
  25         figure(2);fplot(prob,[0,half_prob_heat*5]);
  26         xlabel('heat fraction (1)');ylabel('probability of detection (1)');
  27         title('Logistic fire detection model')
  28         fig
  29
  30         figure(3);fplot(@(T) prob(heat(-T)),[-c,10*c]);
  31         xlabel('time since fire arrival (h)');ylabel('probability of detection (1)');
  32         title('Fire detection probability')
  33         fig
  34
  35         m = 10                  % mesh refinement
  36         dx=300/m, dy=300/m,    % grid step
  37         ix = 10*m, iy=10*m     % grid center indices
  38         nx = 2*ix-1, ny=2*iy-1     % grid dimension
  39         sigma = 2000,     % geolocation error stdev
  40         R = 1            % rate of spread  m/s
  41         R = 3600*R       % rate of spread  m/h
  42
  43         [x,y]=ndgrid([0:nx-1]*dx,[0:ny-1]*dy);  % mesh
  44
  45         ssum = @(a) sum(a(:));    % utility: sum of array
  46
  47         % gaussian weights
  48         d2 = (x(ix,iy)-x).^2 + (y(ix,iy)-y).^2; % distance ^2
  49         w=exp(-d2/(sigma^2));                   % gaussian kernel
  50         w=w/sum(w(:));                          % weights add up to 1
  51
  52         %********************************
  53         tign = @(T)T+(x-x(ix,iy))/R;            % fire arrival time
  54         prbx = @(T)prob(heat(tign(T)));         % detection probability
  55         pgeo = @(T)ssum(w.*prob(heat(tign(T))));% weighted by geolocation error
  56         ageo = @(T)arrayfun(pgeo,T);            % same for array argument
  57         %********************************
  58         figure(4);fplot(@(T) log(ageo(-T)),[-c,10*c]);
  59         xlabel('time since fire arrival (h)');ylabel('log probability of detection (1)');
  60         title('Log likelihood with geolocation error')
  61         fig
  62
  63
  64
  65
  66
  67  %end