-
Notifications
You must be signed in to change notification settings - Fork 5
/
mctaskmeltrejectr.m
64 lines (50 loc) · 1.87 KB
/
mctaskmeltrejectr.m
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
function [simaverage, simerror, simratio] = mctaskmelt(data, prob, uncert, meltdata, Liquid_mass, T, P, binedges, nbins)
samplerows=size(data,1);
sdata=NaN(samplerows,size(data,2));
i=1;
while i<samplerows
% select weighted sample of data
r=rand(length(prob),1);
tempdata=data(prob>r,:);
if i+size(tempdata,1)-1<=samplerows
sdata(i:i+size(tempdata,1)-1,:)=tempdata;
else
sdata(i:end,:)=tempdata(1:samplerows-i+1,:);
end
i=i+size(tempdata,1);
end
% Randomize variables over uncertainty interval
sdata=sdata+sdata.*repmat(uncert,size(sdata,1),1).*randn(size(sdata));
% Randomize ages over uncertainty interval
r=randn(size(sdata(:,1)));
ages=sdata(:,2)+r.*sdata(:,1)/2;
% Reshape composition array
composition=NaN(size(sdata,2)-2,size(sdata,1));
for i=1:(size(sdata,2)-2)
composition(i,:)=sdata(:,i+2)';
end
% Find best-fitting melt percentages
index=zeros(size(composition,2),1);
residual=zeros(size(composition,2),1);
for i=1:size(composition,2)
a=repmat(composition(:,i),1,size(meltdata,2));
[residual(i), index(i)]=min(nansum(((a-meltdata)./a).^2,1));
end
% Take only the half of the samples with the lowest residuals
t = residual<median(residual);
index=index(t);
ages=ages(t);
% Populate temporary variables for each of the elements/ratios of interest
param=NaN(size(index,1),1,3);
param(:,1,1)=Liquid_mass(index); %Meltestimate
param(:,1,2)=T(index); %Tempestimate
param(:,1,3)=potentialtemperature(T(index)+273.15,P)-273.15; %Mantle potential tempestimate
% Find average values for each quantity of interest for each time bin
simaverage=NaN(1,nbins,size(param,3));
simerror=simaverage;
simratio=size(sdata,1)./size(data,1).*sum(t)./length(t);
for j=1:nbins
simaverage(1,j,:)=nanmean(param(ages>binedges(j) & ages<binedges(j+1),:,:),1);
simerror(1,j,:)=nansem(param(ages>binedges(j) & ages<binedges(j+1),:,:),1);
end
end