mcmc算法实现代码

Question

你好你好，我刚刚接触代码，请问这样子的题目怎么在matlab中用代码实现呀？
数据如下：

0 0.2
0.08 17.5
0.17 31.2
0.25 45.3
0.33 62.7
0.42 78.6
0.5 103
0.58 123.5
0.67 133.2
0.75 171.9
0.83 178.1
0.92 164.9
1 174.3
1.08 197.4
1.17 220.2
1.25 231.6
1.33 265.7
1.42 279
1.5 300.7
1.58 299
1.67 294
1.75 278.1
1.83 363.8
1.92 365.7
2 409
2.08 408.3
2.17 359.5
2.25 438.3
2.33 413.6
2.42 459.5
2.5 428.3
2.58 465.7
2.67 429.3
2.75 454.6
2.83 498.6
2.92 449.1
3 497.4
3.08 496.2
3.17 507.9
3.25 466.1
3.33 665.8
3.42 638
3.5 576.9
3.58 581.2
3.67 538.6
3.75 620.5
3.83 624
3.92 561.1
4 595.2
4.08 570.4
4.17 627.5
4.25 704
4.33 544.4
4.42 670.1
4.5 606.6
4.58 622.8
4.67 673.1
4.75 720.3
4.83 742
4.92 664.7
5 794
5.08 722.5
5.17 789.4
5.25 751.3
5.33 634.4
5.42 792.5
5.5 792.7
5.58 839.6
5.67 800
5.75 661.5
5.83 903
5.92 812.3
6 821.5
6.08 784.1
6.17 828
6.25 785.6
6.33 877.8
6.42 903.1
6.5 839
6.58 917.5
6.67 756.1
6.75 896.6
6.83 871.7
6.92 966.7
7 923.2
7.08 748.6
7.17 855.1
7.25 937.7
7.33 887.2
7.42 942.9
7.5 840.5
7.58 859.7
7.67 952.2
7.75 1030.9
7.83 1017.1
7.92 959.5
8 927.5
8.08 1054
8.17 831.4
8.25 800.9
8.33 1023.2
8.42 817.3
8.5 1110.5
8.58 1072.6
8.67 1042.4
8.75 967.3
8.83 896.9
8.92 1060
9 974.1
9.08 1020
9.17 859.5
9.25 945.9
9.33 903.3
9.42 889.7
9.5 1086.7
9.58 886.9
9.67 972.5
9.75 1015.5
9.83 850.3
9.92 981.9
10 1196
10.08 1092.4
10.17 1067.4
10.25 1027.4
10.33 1054.4
10.42 1026.8
10.5 1042.7
10.58 981.4
10.67 927.2
10.75 985.7
10.83 1151
10.92 1119.4
11 1182.3
11.08 1174.8
11.17 932.5
11.25 1063.5
11.33 809.7
11.42 1035.5
11.5 919.9
11.58 1058.6
11.67 1104
11.75 971.1
11.83 1052.8
11.92 974.2
12 1038.1
12.08 1092.5
12.17 1032.5
12.25 895
12.33 931.8
12.42 1044.6
12.5 1067.1
12.58 1080.7
12.67 1117.3
12.75 914.9
12.83 907
12.92 1039.6
13 1120.4
13.08 1201.9
13.17 1149.3
13.25 988
13.33 1073.9
13.42 933.9
13.5 1336
13.58 1029
13.67 997.4
13.75 1058.5
13.83 1080.5
13.92 953.1
14 1169.6
14.08 1033.7
14.17 1043.7
14.25 1155.6
14.33 976.1
14.42 1139.8
14.5 845.2
14.58 884.3
14.67 1147
14.75 1052
14.83 1211.4
14.92 1172.2
15 1042.4
15.08 1060.2
15.17 1191
15.25 1217.2
15.33 1242.9
15.42 1189.6
15.5 1189.2
15.58 1104.9
15.67 1233.6
15.75 1180.5
15.83 1187.3
15.92 1132.1
16 1222.7
16.08 1167
16.17 1271
16.25 1014.9
16.33 1053.8
16.42 1167.4
16.5 1161.1
16.58 1125.7
16.67 1242
16.75 1219.3
16.83 913.4
16.92 860.8
17 1056.9
17.08 939.4
17.17 1097
17.25 1057.8
17.33 964.9
17.42 908
17.5 921.4
17.58 1138.6
17.67 1072.4
17.75 1035.5
17.83 1120.6
17.92 893.2
18 989.9
18.08 1005.2
18.17 1147.3
18.25 1028.9
18.33 1248.8
18.42 911
18.5 1053.2
18.58 1167.9
18.67 1238.8
18.75 1021.9
18.83 1041.5
18.92 1045.4
19 1009.3
19.08 1029.9
19.17 912.5
19.25 1047.9
19.33 1049.2
19.42 964.5
19.5 1131.9
19.58 1145.2
19.67 850.5
19.75 1282.9
19.83 1017
19.92 1093
20 912.6
20.08 1253.8
20.17 1129.8
20.25 1087.7
20.33 945.7
20.42 1195.6
20.5 1055.6
20.58 1073.5
20.67 1146.8
20.75 1285.6
20.83 948.1
20.92 1094.2
21 1138.9
21.08 952.6
21.17 944.1
21.25 1068.5
21.33 708.9
21.42 953.8
21.5 880.2
21.58 1021
21.67 1056.2
21.75 1020.7
21.83 1070.5
21.92 1130.4
22 1000.4
22.08 874.8
22.17 1032
22.25 1069.8
22.33 974
22.42 1049.5
22.5 1108.4
22.58 1008.9
22.67 1037.1
22.75 911.8
22.83 983.2
22.92 1189
23 1077.2
23.08 918.4
23.17 1040
23.25 1113.9
23.33 1098.2
23.42 798.3
23.5 883.9
23.58 880.9
23.67 966.6
23.75 972.4
23.83 1154.4
23.92 1004.1
24 1053.8
24.08 1104.7
24.17 1234.2
24.25 1161.9
24.33 997.4
24.42 941.7
24.5 920.6
24.58 1038.5
24.67 1050.2
24.75 887.5
24.83 884.2
24.92 853.5
25 1063.2

这个是我根据助教修改的代码，但是不明白每一区域的原理，请问怎么修改下面的代码呀？

GFPData = csvread("D:\上科大\高级统计学\第四次作业\transfection.csv");

% % % Measurement errors for the calculation of likelihhood function
errorGFP = 50; % % measurement error for GFP

% % % Set up initial conditions
NIter = 301;
t = 1;
theta0 = [500,0.2,1,50,30]; % % [m0,G0,k,a,b]
thetaCurrent = theta0;

priorG = [100,2];

x = NaN(NIter,2); % % % sampled theta
L = NaN(NIter,1); % % % Posterior P (prior + likelihood);
A = false(NIter,1); % % % 1 for accept

x(1,:) = theta0;

logLikeCurrent = logLike(thetaCurrent,GFPData,errorGFP);
logPriorCurrent = logPrior(thetaCurrent,priorG);
L(1,:) = logPriorCurrent + logLikeCurrent;

% % % Use mvnrnd(log(thetaCurrent),S) to propose log(thetaNext)
toLog = @(theta) log([-log(theta(1)),theta(2:301)]);
backLog =@(thetaLog) exp([-exp(thetaLog(1)),thetaLog(2:301)]);

meanAcceptRate = 0.3;
S_scaling = 1;
S = eye(length(theta0)); % % % covariance for theta
M = toLog(theta0);

while t<NIter
    t = t+1;
    
    % % % Propose thetaNextLog according to thetaCurrentLog
    thetaNext = backLog(mvnrnd(toLog(thetaCurrent),S*S_scaling,1));

    % % % Calculate likelihood and priors for thetaNext
    logLikeNext = logLike(thetaNext,GFPData,errorGFP);
    logPriorNext = logPrior(thetaNext,priorG);

    % % % Get the acceptance rate
    acceptRate = logPriorNext + logLikeNext - logPriorCurrent - logLikeCurrent ;
    acceptRate = min(1,exp(acceptRate));

    % % % Accept with probability=acceptRate
    if rand<acceptRate % % accept the proposed theta
        thetaCurrent = thetaNext;
        x(t,:) = thetaNext;
        logPriorCurrent = logPriorNext;
        logLikeCurrent = logLikeNext;
        L(t,1) = logPriorNext + logLikeNext;
        A(t) = true;
    else % % reject the proposed theta and accept current theta
        x(t,:) = thetaCurrent;
        L(t,1) = logPriorCurrent + logLikeCurrent;
    end

    % % % Update covariance
    y = 1/3/sqrt(t);
    S_scaling = exp(log(S_scaling) + y*(acceptRate - meanAcceptRate)); % % % scaling factor for sigma
    S = S + y*((toLog(thetaCurrent) - M)'*(toLog(thetaCurrent) - M) - S);
    M = M + y*(toLog(thetaCurrent) - M);

end


function logLikelihood = logLike(theta,GFPData,errorGFP)
% % % 根据作业题修改似然函数，由于作业中只测量了GFP的荧光强度，并不知道mRNA的测量值，所以只用根据当前的模型参数计算GFP强度的似然
    [~,AU] = ode45(@(t,au) tCourse(t,au,theta(3),theta(4)), [0,GFPData(:,1)'], [theta(1);theta(2)]);
    logLikelihood =  ((AU(2:end,1)-GFPData(:,2))/errorGFP).^2;
    logLikelihood = -sum(logLikelihood(:))/2;
end


function prior = ((theta(2)-priorG(1))/priorG(2)).^2 + ...   
                 ((theta(3)-priorAlpha(1))/priorAlpha(2)).^2 + ...
                 ((theta(4)-priorBeta(1))/priorBeta(2)).^2;
prior = -prior/2;
priorAlpha = 1/theta(3); % alpha的无信息先验
priorBeta = 1/theta(4); % beta的无信息先验
prior = -((theta(2)-priorG(1))/priorG(2)).^2/2 - log(priorAlpha) - log(priorBeta); 

end

function dGdt = tCourse(~, measure,alpha,beta)

% % % 根据作业题修改方程
     % % % measure(1) is M and measure(2) is Glc
         dGdt = [  alpha*(measure(1)/(1+measure(1)))*measure(1);...
               -alpha/beta*(measure(2)/(1+measure(2)))*measure(1)];
end

Answer 1

基于Monster 组和AIGC的调写：若有帮助，还望采纳~~：

修改后：

GFPData = readtable('C:\Users\Monster\Desktop\新建文件夹\matlab\transfection.csv', 'ReadVariableNames', false);
GFPData = table2array(GFPData);


 
% % % Measurement errors for the calculation of likelihhood function
errorGFP = 50; % % measurement error for GFP
 
% % % Set up initial conditions
NIter = 301;
t = 1;
theta0 = [500,0.2,1,50,30]; % % [m0,G0,k,a,b]
thetaCurrent = theta0;
 
priorG = [100,2];
 
x = NaN(NIter, length(theta0));  % 这将创建一个有 NIter 行和与 theta0 长度相同列数的矩阵
x(1,:) = theta0;  

L = NaN(NIter,1); % % % Posterior P (prior + likelihood);
A = false(NIter,1); % % % 1 for accept
 
x(1,:) = theta0;
 
logLikeCurrent = logLike(thetaCurrent,GFPData,errorGFP);
logPriorCurrent = logPrior(thetaCurrent,priorG);
L(1,:) = logPriorCurrent + logLikeCurrent;
 
% % % Use mvnrnd(log(thetaCurrent),S) to propose log(thetaNext)
toLog = @(theta) log(theta);

backLog = @(thetaLog) exp(thetaLog);

 
meanAcceptRate = 0.3;
S_scaling = 1;
S = eye(length(theta0)); % % % covariance for theta
M = toLog(theta0);
 
while t<NIter
    t = t+1;
    
    % % % Propose thetaNextLog according to thetaCurrentLog
    thetaNext = backLog(mvnrnd(toLog(thetaCurrent),S*S_scaling,1));
 
    % % % Calculate likelihood and priors for thetaNext
    logLikeNext = logLike(thetaNext,GFPData,errorGFP);
    logPriorNext = logPrior(thetaNext,priorG);
 
    % % % Get the acceptance rate
    acceptRate = logPriorNext + logLikeNext - logPriorCurrent - logLikeCurrent ;
    acceptRate = min(1,exp(acceptRate));
 
    % % % Accept with probability=acceptRate
    if rand<acceptRate % % accept the proposed theta
        thetaCurrent = thetaNext;
        x(t,:) = thetaNext;
        logPriorCurrent = logPriorNext;
        logLikeCurrent = logLikeNext;
        L(t,1) = logPriorNext + logLikeNext;
        A(t) = true;
    else % % reject the proposed theta and accept current theta
        x(t,:) = thetaCurrent;
        L(t,1) = logPriorCurrent + logLikeCurrent;
        disp(['Iteration ', num2str(t), ': theta = ', num2str(thetaCurrent)]);

    end
    % 绘制参数 k 的后验分布
    histogram(x(:,3)); % 假设 k 是 theta 中的第三个参数
    title('Posterior Distribution of Parameter k');
 
    % % % Update covariance
    y = 1/3/sqrt(t);
    S_scaling = exp(log(S_scaling) + y*(acceptRate - meanAcceptRate)); % % % scaling factor for sigma
    S = S + y*((toLog(thetaCurrent) - M)'*(toLog(thetaCurrent) - M) - S);
    M = M + y*(toLog(thetaCurrent) - M);
 
end

figure;
for i = 1:4
    subplot(2, 2, i);
    plot(x(:,i));
    xlabel('Step');
    ylabel(['Parameter ', num2str(i)]);
    title(['Parameter ', num2str(i), ' vs Steps']);
end
% 绘制直方图和拟合曲线
figure;
for i = 1:4
    subplot(2, 2, i);
    h = histogram(x(:,i));
    hold on;
    % 检查数据
    binCenters = h.BinEdges(1:end-1) + diff(h.BinEdges) / 2;
    counts = h.Values;

    % 确保数据有效
    validIndices = ~isinf(counts) & ~isnan(counts);
    binCenters = binCenters(validIndices);
    counts = counts(validIndices);

    % 尝试拟合一次多项式模型
    if ~isempty(binCenters) && ~isempty(counts)
        fitResults = fit(binCenters.', counts.', 'poly1');
        plot(fitResults, binCenters, counts);
    end
    xlabel(['Parameter ', num2str(i)]);
    ylabel('Count');
    title(['Histogram and Fit for Parameter ', num2str(i)]);
    hold off;
end

 
function logLikelihood = logLike(theta,GFPData,errorGFP)
    [~,G] = ode45(@(t,G) tCourse(t,G,theta(3),theta(4),theta(5)), GFPData(:,1), [theta(1);theta(2)]);
    logLikelihood = -sum(((G(:,2)-GFPData(:,2))/errorGFP).^2)/2;
end
 
 
function prior = logPrior(theta,priorG,priorAlpha,priorBeta)
    prior = -((theta(2)-priorG(1))/priorG(2)).^2/2 - log(1/theta(3)) - log(1/theta(4)) - log(1/theta(5)); 
end

 
function dGdt = tCourse(t, G, k, alpha, beta)
    m = G(1);
    GFP = G(2);
    dm = -alpha * m;
    dGFP = k * m - beta * GFP;
    dGdt = [dm; dGFP];
end