DevelopmentalExpressionModeling/GetWithErrorBars.m at master · DynamicsAndNeuralSystems/DevelopmentalExpressionModeling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
function GetWithErrorBars(numRepeats)
% Lots of simulations of autocorrelated spatial maps
%-------------------------------------------------------------------------------

if nargin < 1
    numRepeats = 50;
end
%-------------------------------------------------------------------------------

params = GiveMeDefaultParams();

numTimePoints = 7;
lambdaEst = zeros(numRepeats,numTimePoints);
lambdaErr = zeros(numRepeats,numTimePoints);
strengthEst = zeros(numRepeats,numTimePoints);
strengthErr = zeros(numRepeats,numTimePoints);
offSetEst = zeros(numRepeats,numTimePoints);
offSetErr = zeros(numRepeats,numTimePoints);

for i = 1:numRepeats
    fprintf(1,'~~~~~~~ITERATION %u/%u~~~~~~~\n',i,numRepeats);
    [maxExtent,paramStruct] = simulateGrid(params,false);
    paramNames = coeffnames(paramStruct{1});

    % Mean estimates:
    lambdaEst(i,:) = cellfun(@(x)1/x.n,paramStruct);
    strengthEst(i,:) = cellfun(@(x)x.A,paramStruct);
    offSetEst(i,:) = cellfun(@(x)x.B,paramStruct);

    % CIs of estimates:
    CIs = cellfun(@(x)confint(x),paramStruct,'UniformOutput',false);
    % Invert for lambda:
    nIndex = strcmp(paramNames,'n');
    for j = 1:numTimePoints
        CIs{j}(:,nIndex) = 1./CIs{j}(:,nIndex);
    end
    % Store errors:
    strengthErr(i,:) = cellfun(@(x)diff(x(:,1))/2,CIs);
    offSetErr(i,:) = cellfun(@(x)diff(x(:,2))/2,CIs);
    lambdaErr(i,:) = cellfun(@(x)diff(x(:,3))/2,CIs);
end

%-------------------------------------------------------------------------------
% Just in case
save('HailMary.mat')

%-------------------------------------------------------------------------------
f = figure('color','w');
myColors = BF_getcmap('dark2',7,1);
for i = 1:3
    subplot(1,3,i); hold('on')
    switch i
    case 1
        paramEst = lambdaEst;
        paramErr = lambdaErr;
        ylabel('Spatial scale, \lambda')
    case 2
        paramEst = strengthEst;
        paramErr = strengthErr;
        ylabel('Strength, A')
    case 3
        paramEst = offSetEst;
        paramErr = offSetErr;
        ylabel('Offset, f_0')
    end

    % Fit linear regression:
    if i==1
        ft = fittype('c*x');
        [c,Stats] = fit(maxExtent,mean(paramEst)',ft);
        f_handle = @(x) c.c*x;
        xRange = linspace(0,maxExtent(end));
        plot(xRange,f_handle(xRange),':k')
    end

    %-------------------------------------------------------------------------------
    % Plot (colored) parameter variation
    params = struct();
    params.doSubsample = true;
    params.thisBrainDiv = 'brain';
    params.thisCellType = 'allCellTypes';
    params.includeAdult = false;

    [params,fittedParams,CIs,goodTimePoint] = LoadParameterFits(params);

    lambdaEmpirical = cellfun(@(x)1/x.n,fittedParams);
    lambdaErrs = cellfun(@(x)diff(x(:,3))/2,CIs);
    strengthEmpirical = cellfun(@(x)x.A,fittedParams);
    strengthErrs = cellfun(@(x)diff(x(:,1))/2,CIs);
    f0Empirical = cellfun(@(x)x.B,fittedParams);
    f0Errs = cellfun(@(x)abs(diff(x(:,2)))/2,CIs);
    paramMeanValues = [lambdaEmpirical,strengthEmpirical,f0Empirical];
    paramErrValues = [lambdaErrs,strengthErrs,f0Errs];

    for t = 1:numTimePoints
        % Simulation:
        errorbar(maxExtent(t),mean(paramEst(:,t)),mean(paramErr(:,t)),'o','color',...
                                myColors{t},'LineWidth',1.8)
        % Empirical data:
        empirical = paramMeanValues(:,i);
        empiricalErrs = paramErrValues(:,i);
        smallOffset = 0.1;
        errorbar(maxExtent(t)+smallOffset,empirical(t),empiricalErrs(t),'o',...
                            'color',brighten(myColors{t},+0.5),'LineWidth',1.8)
    end

    xlabel('Brain size, dmax (mm)')
end
f.Position = [1000        1116         831         222];

%-------------------------------------------------------------------------------
% Save to file:
fileName = sprintf('Model_Simo_%uReps.svg',numRepeats);
saveas(f,fileName,'svg')
fprintf(1,'Saved to %s\n',fileName);

end