diff --git a/run.m b/run.m
index a90264e..b79e315 100644
--- a/run.m
+++ b/run.m
@@ -302,4 +302,48 @@ phaseLoadQError=reshape(phaseLoadQError,size(phaseLoadQError,1)*size(phaseLoadQE
 % figure()
 % barPlot( phaseLoadPError, 10,['相对误差%'],['分布密度'],['有功负荷误差']);
 % figure()
-% barPlot( phaseLoadQError, 10,'相对误差%','分布密度','无功负荷误差');
\ No newline at end of file
+% barPlot( phaseLoadQError, 10,'相对误差%','分布密度','无功负荷误差');
+%% 计算统计量
+%目标函数均值
+rVoltABCV=conj([VoltpA';VoltpC';VoltpC']);
+mVoltABCV=conj([mVoltpA';mVoltpB';mVoltpC']);
+rPD3P=real([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
+    phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
+    phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
+    ]);
+rPD3P=rPD3P(:,Loadi);
+mPD3P=real([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
+    mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
+    mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
+    ]);
+mPD3P=mPD3P(:,Loadi);
+PD3P=real(SEThreeLoad);
+rQD3P=imag([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
+    phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
+    phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
+    ]);
+rQD3P=rQD3P(:,Loadi);
+mQD3P=imag([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
+    mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
+    mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
+    ]);
+mQD3P=mQD3P(:,Loadi);
+QD3P=imag(SEThreeLoad);
+JMeasurement=sum(sum((( abs(mVoltABCV)-abs(SEVoltpABC) )./ abs(mVoltABCV)./sigma).^2))+sum(sum(((mPD3P-PD3P)./mPD3P./sigma).^2))+sum(sum(((mQD3P-QD3P)./mQD3P./sigma).^2));
+%估计误差统计
+%量测量数量
+Busnum=busNum;
+mCount=Busnum*3+length(Loadi)*3*2;
+%估计量质量
+%有PD，QD的
+SEMeasurement=sum(sum((( abs(rVoltABCV)-abs(SEVoltpABC) )./abs(mVoltABCV)./sigma).^2))+sum(sum(((rPD3P-PD3P)./mPD3P./sigma).^2))+sum(sum(((rQD3P-QD3P)./mQD3P./sigma).^2));
+SEMeasurement=(SEMeasurement/mCount)^.5;
+%没有PD，QD的
+SEMeasurementNOPQ=sum(sum(((abs(rVoltABCV)-abs(SEVoltpABC))./abs(mVoltABCV)./sigma).^2));
+SEMeasurementNOPQ=(SEMeasurementNOPQ/mCount)^.5;
+%量测量质量
+SEMeasurementM=sum(sum((( abs(rVoltABCV)-abs(mVoltABCV) )./abs(mVoltABCV)./sigma).^2))+sum(sum(((rPD3P-mPD3P)./mPD3P./sigma).^2))+sum(sum(((rQD3P-mQD3P)./mQD3P./sigma).^2));
+SEMeasurementM=(SEMeasurementM/mCount)^.5;
+%没有PD，QD的
+SEMeasurementNOPQM=sum(sum((( abs(rVoltABCV)-abs(mVoltABCV))./abs(mVoltABCV)./sigma).^2));
+SEMeasurementNOPQM=(SEMeasurementNOPQM/mCount)^.5;
\ No newline at end of file