机器算法验证 - 如何比较不同测量仪器的测量值和不确定度？ - 吾爱随笔录

如何比较不同测量仪器的测量值和不确定度？

机器算法验证 r 分布重复测量测量误差测量

2022-04-01 02:57:16

我有两个不同的测量仪器，A 和 B，都测量相同的物理特性 $x$ 一个物体但具有不同的“质量”：B给出了一个已知不确定度的测量值，而我不知道A给出的测量值的不确定度。

我有 $N$ 不同的物体，我测量属性 $x$ 对于所有带A的人，所以我得到了一个测量列表 $L_A=\{x_{A1}, x_{A2},\ldots,x_{AN}\}$ 在哪里 $x_{Ai}$ 是属性的度量 $x$ 为了 $i$ -th 对象。

对象没有标记，给定一个对象，我只知道它的测量属于 $L_A$ 但我无法从中提取测量值 $L_A$ . 此外，我不能用 A 来测量我过去用 A 测量过的对象。

然后我随机选择 $M \lt N$ 从对象 $N$ 对象。我测量所有 $M$ 带有 B 的样本对象，我得到一个测量列表 $L_B=\{x_{B1}, x_{B2},\ldots,x_{BM}\}$ . 请注意，下标中的索引不是对象的标签，所以我无法直接比较 $x_{A1}$ 和 $x_{B1}$ .

是否可以用上述数据估计 A 给出的测量值的不确定性？

我在考虑比较经验累积分布函数 $L_A$ 与其中之一 $L_B$ 但这只是一个想法，我无法进一步详细说明。

是否有任何既定标准可以涵盖我的问题？例如，我找到了对“ISO 5725：测量方法和结果的准确性（真实性和精确度）”的引用，但我无法访问它。

更新：

我发现我的问题类似于如何测试从两个设备读取的数据是否显着不同？在我读到Michael Lew的答案的地方，他推荐了LUDBROOK, John 的论文。比较测量者和测量方法的统计技术：批判性评论。临床和实验药理学和生理学，2002，29.7：527-536。但不幸的是，在我看来，这篇论文需要测量之间的配对。

更新 2：

我写了一个 R 脚本来模拟我的问题。

set.seed(42)

instrument_measurement <- function(true_value,gain,offset,dispersion)
# The instrument has three parameters: the true_value is transformed by means
# of a linear transformation described by parameters offset and gain; then there
# is a dispersion parameter. An ideal instrument would have gain=1, offset=0
# and dispersion that approaches to zero.
{
  return(rnorm(length(true_value),mean=gain*true_value-offset,sd=dispersion))
}

N=1000
true_mean = 0
true_sd = 1
# I simulate the property to be measured for N objects, here the property is
# normally distributed...
true_values = rnorm(N,mean=true_mean,sd=true_sd)

# but it could also be a mixture of normal distributions:

# components <- sample(1:3,prob=c(1/7,5/7,1/7),size=N,replace=TRUE)
# mus <- c(-0.2,0,+0.3)
# sds <- sqrt(c(0.05,0.05,0.05))
# true_values <- rnorm(n=N,mean=mus[components],sd=sds[components])
# plot(density(true_values))


# The "quality" of instrument B is "good enough" to measure the true values:
gain_B = 1
offset_B = true_sd/10
dispersion_B = true_sd/10

# The instrument B has a lower "quality" than the one of instrument A:
gain_A = 1.1*gain_B
offset_A=-2*offset_B
dispersion_A=5*dispersion_B

# I simulate the measuremente made by instrument A:
L_A = instrument_measurement(true_values,gain_A,offset_A,dispersion_A)

# I make the sample:
sample_to_measure_with_B = sample(true_values,100,replace=F)

# I simulate the measuremente made by instrument B:
L_B = instrument_measurement(sample_to_measure_with_B,gain_B,offset_B,dispersion_B)

# I plot the empirical CDF of the true values, of the measurements made with
# instrument A and of the measurements made with instrument B
plot(ecdf(true_values),col="grey",main="",xlab="x, measured property",ylab="value of empirical CDF")
lines(ecdf(L_A),col="blue")
lines(ecdf(L_B),col="orange")
legend(x=(max(true_values)+mean(true_values))/2,y=.5,legend=c("true","A","B"),col=c("grey","blue","orange"),lty=c(1,1,1))
title("Empirical CDFs")

参考脚本，可以估计gain_A,offset_A和dispersion_AfromL_A和L_B? 估计中的不确定性是什么？

我有一个定义成本函数并尝试在参数空间中最小化它的不优雅的想法gain，offset并且dispersion：

function <- ecdf_distance(ecdf1,ecdf2)
{
# return 0 if ecdf1 is "equal" to ecdf2
# return a positive scalar that measures the difference between ecdf1 and ecdf2
}

function <- cost(parameters)
{
L = instrument_measurement(L_B,parameters$gain,parameter$offset,parameter$dispersion)
return(ecdf_distance(ecdf(L_A),ecdf(L))
}

我做了一些测试gain=1，但没有运气......成本函数似乎是恒定的dispersion......我担心我缺乏一些关于这个问题的理论/数学:-)

在此处输入图像描述

2个回答

我更像是一个工程人员，而不是统计人员，所以让我们谈谈具体细节。

寻找：

估计的对象度量不确定性（状态协方差）
（可能）用于发现测量不确定度的既定标准测量。

鉴于：

每次测量后随机更换被测对象，不允许配对比较。
具有不同测量不确定度的两种测量工具，其中一种不确定度没有被表征。

注意事项：

校准。大多数测量系统都根据高质量标准进行校准，因此它们的偏差最小。我应该能够假设这一点，但是您的系统表明偏见是一个问题。
样本量。进行 100 次测量非常弱。我更喜欢至少 300，但这与我的数据有关。您需要确保获得足够的样本以最大限度地减少估计中的误差，但不要太多以至于您被埋没。我可以轻松获得 1000 万个样本，但它超出了 MatLab 和我的笔记本电脑的范围，除了基础知识之外还有那么多行。

方法：

将您的代码从 R 转换为我熟悉的语言 MatLab。
查看 CDF 域插值和散点图。

您的代码的 MatLab：

function MySimulation

%houekeeping
clc;

%parameters
N=1000
true_mean = 0
true_sd = 1
% I simulate the property to be measured for N objects, here the property is
% normally distributed...
true_values = normrnd(true_mean,true_sd,N,1);

% but it could also be a mixture of normal distributions:
% it could be anything.  A Gaussian mixture is crazy-tame compared to what
% it could be, but we have to assume something to start this.

% The "quality" of instrument B is "good enough" to measure the true values:
gain_B = 1;
offset_B = true_sd/10;
dispersion_B = true_sd/10;

% The instrument B has a lower "quality" than the one of instrument A:
gain_A = 1.1*gain_B;
offset_A=-2*offset_B;
dispersion_A=5*dispersion_B;

% I simulate the measuremente made by instrument A:
L_A = instrument_measurement(true_values,gain_A,offset_A,dispersion_A);

% I make the sample:
% sample_to_measure_with_B = sample(true_values,100,replace=F)
sample_to_measure_with_B = randsample(true_values,100);

% I simulate the measuremente made by instrument B:
L_B = instrument_measurement(sample_to_measure_with_B,gain_B,offset_B,dispersion_B);

% I plot the empirical CDF of the true values, of the measurements made with
% instrument A and of the measurements made with instrument B

figure(1); clf; hold on
h1=cdfplot(true_values);
set(h1,'Color',0.5.*[1,1,1],'Linewidth',2);

h2=cdfplot(L_A);
set(h2,'Color','b','Marker','+');

h3=cdfplot(L_B);
set(h3,'Color',[1,0.65,0],'Marker','o');

legend('true','A','B','location','northwest')
title('Empirical CDFs')
xlabel('x-measured property')
ylabel('value of empirical CDF')


function [out] = instrument_measurement(true_value,gain,offset,dispersion)
% The instrument has three parameters: the true_value is transformed by means
% of a linear transformation described by parameters offset and gain; then there
% is a dispersion parameter. An ideal instrument would have gain=1, offset=0
% and dispersion that approaches to zero.
out=normrnd(gain*true_value-offset,dispersion,length(true_value),1);
return

你的图片，翻译：在此处输入图像描述

10 倍更高的采样率：在此处输入图像描述

现在，当我们使用以下方法绘制 A 与 B 时

[xa,fa]=ecdf(L_A);
[xb,fb]=ecdf(L_B);

xb2=xa;
fb2=interp1(xb,fb,xb2);

我们得到这个：

在此处输入图像描述

如果你看它，你会发现，除了尾部失控之外，这两个分布之间的关系基本上是线性的。它是比例和偏移量。我会使用样本量来确定在哪里截断尾部，然后通过它拟合一条分析线。然后我可以将一个传感器读数转换为另一个。

通过 CDF 映射的技巧是如何处理许多分布。您需要确保您的样本量足够大以捕捉潜在分布的特征。

如果您一次处理一口，而不是一次处理 10k 个样本，那么您可以使用卡尔曼滤波器来确定均值和方差，或 A 和 B 之间的变换。前面是“rev1”蛮力强制的方法。这是不雅的。它有漏洞和弱点。它没有强烈的特征。它没有建立在惊人的强大理论基础上。然而，它相当快且“足够好”，这也是衡量善良的好方法。

关于韦尔奇有一个很好的介绍。（链接）

一个关键词是“传感器融合”。这是有关该主题的第二个链接。（链接）

祝你好运。

您用来“模拟您的问题”的模型几乎可以逐字地使用贝叶斯估计来估计您感兴趣的参数。这是我将使用的模型（使用与您相同的符号）：

L_{B} \sim N o r m a l (μ, σ) x_{i} \sim N o r m a l (μ, σ) f o r i f r o m 1 t o N L_{A i} \sim N o r m a l (x_{i} \cdot g a i n - o f f s e t, d i s p e r s i o n) f o r i f r o m 1 t o N

$L_B \sim \mathrm{Normal}(\mu, \sigma) \\ x_i \sim \mathrm{Normal}(\mu, \sigma) \mathrm{\ for\ i\ from\ 1\ to\ N} \\ L_{Ai} \sim \mathrm{Normal}(x_i \cdot \mathrm{gain} - \mathrm{offset}, \mathrm{dispersion}) \mathrm{\ for\ i\ from\ 1\ to\ N} \\$

与您的问题相比，此模型中明显的遗漏是我不包括假设某些相同 $x_i$ 由 B 测量的 s 然后可以由 A 再次测量。这可能会被添加，但我不完全确定如何添加。

该模型在下面的 R & JAGS 中使用非常模糊、几乎平坦的先验实现，使用的数据是您在问题中生成的数据：

library(rjags)

model_string <- "model{
  for(i in 1:length(L_B)) {
  L_B[i] ~ dnorm(mu, inv_sigma2) # <- reparameterizing sigma into precision 
                                 #    needed because of JAGS/BUGS legacy.  
  }
  for(i in 1:length(L_A)) {
    x[i] ~ dnorm(mu, inv_sigma2)
    L_A[i] ~ dnorm(gain * x[i] - offset , inv_dispersion2)
  }

  mu ~ dnorm(0, 0.00001)
  inv_sigma2 ~ dgamma(0.0001, 0.0001) 
  sigma <- sqrt(1 / inv_sigma2)
  gain ~ dnorm(0, 0.00001) T(0,)
  offset ~ dnorm(0, 0.00001)
  inv_dispersion2 ~ dgamma(0.0001, 0.0001)
  dispersion <- sqrt(1 / inv_dispersion2)
}"

让我们运行它，看看它的效果如何：

model <- jags.model(textConnection(model_string), list(L_A = L_A, L_B = L_B), n.chains=3)
update(model, 3000)
mcmc_samples <- coda.samples(model, c("mu", "sigma", "gain", "offset", "dispersion"), 200000, thin=100)
apply(as.matrix(mcmc_samples), 2, quantile, c(0.025, 0.5, 0.975))
##       dispersion   gain      mu   offset  sigma
## 2.5%     0.01057 0.1366 -0.3116 -0.51836 0.9365
## 50%      0.18657 1.0745 -0.1099 -0.26950 1.0675
## 97.5%    1.20153 1.2846  0.1051 -0.04409 1.2433

结果估计值与您在生成数据时使用的值相当接近：

c(gain_A, offset_A, dispersion_A)
## [1]  1.1 -0.2  0.5

...除了，也许，分散。但是随着更多的数据，可能更知情的先验和运行 MCMC 采样的时间更长，这个估计应该更好。

其它你可能感兴趣的问题

上一篇选择要包含在回归中的稀疏主成分的数量下一篇关于将具有不同反应量表的项目组合成复合量表的建议/文献？