利用python或Matlab对非正态分布的归一化时间序列矩阵进行对数转换、正态转换、正交转换为独立的标准正态多变量,对蒙特卡洛模拟的正态分布结果进行逆转换

#利用python或Matlab对非正态分布的归一化时间序列矩阵进行对数转换、正态转换、正交转换为独立的标准正态多变量,利用上述矩阵转换过程对具有相同时长的正态分布的时间序列进行逆转换。
初始非正态分布时间序列矩阵类似下图:

img


数据获取为以下链接
链接:https://pan.baidu.com/s/1_-f1EJn0raHb_apJK82wzA?pwd=fcno
提取码:fcno

具体的处理过程如下图所示:

img

img

img

各处理过程结果的示例参考上一篇问题
https://ask.csdn.net/questions/7960595?spm=1001.2014.3001.5501

img


matlab可能会更好的解决这个问题

在Python中,使用NumPy和SciPy库来进行对数转换、正态转换和正交转换。下面是一个小小示例代码:

import numpy as np  
from scipy.stats import multivariate_normal  
  
# 生成一个非正态分布的随机矩阵  
np.random.seed(0)  
X = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], size=100)  
  
# 对数转换  
logX = np.log(X)  
  
# 正态转换  
X_norm = np.exp(logX - np.tile(np.mean(logX, axis=0), (len(logX), 1)))  
  
# 正交转换  
X_ort = np.dot(X_norm, np.linalg.inv(np.cov(X_norm)))  
  
# 蒙特卡洛模拟的正态分布结果  
mu = np.array([1, 2])  
cov = np.array([[1, 0.5], [0.5, 1]])  
Y = multivariate_normal.rvs(mean=mu, cov=cov, size=100)  
  
# 逆转换  
Y_ort = np.dot(Y, np.linalg.inv(np.cov(Y)))  
Y_norm = np.exp(np.dot(Y_ort, np.dot(np.linalg.inv(np.cov(X_norm)), np.tile(np.mean(logX, axis=0), (len(Y_ort), 1)))) - np.tile(np.mean(logX, axis=0), (len(Y_ort), 1)))  
Y_log = np.exp(np.dot(Y_norm, np.dot(np.linalg.inv(np.cov(logX)), np.tile(np.mean(logX, axis=0), (len(Y_norm), 1)))) - np.tile(np.mean(logX, axis=0), (len(Y_norm), 1)))

#未完待续,如有帮助,恭请采纳

基于ChatGPT4与博主叶秋学长的回答,望采纳!!!有其他问题也可以询问我哦💕:

首先,你需要对数据进行对数转换。这可以通过numpy库的log函数来实现:

import numpy as np

# 假设 data 是你的数据矩阵
log_data = np.log(data)

然后,为了使数据正态化,你可以使用scipy库的boxcox函数。这个函数可以找到一个合适的lambda值,将数据转换为正态分布:

from scipy import stats

# 注意,boxcox函数只接受正值。如果你的数据中有负值或零,你需要先将它们转换为正值。
data_positive = log_data - np.min(log_data) + 1
normalized_data, best_lambda = stats.boxcox(data_positive)

接下来,你需要将数据进行正交化。这可以通过numpy的linalg.qr函数来实现,它可以进行QR分解,得到一个正交矩阵:

Q, R = np.linalg.qr(normalized_data)

这样,你得到的Q就是一个正交矩阵,它的每一列都是独立的标准正态分布的。

然后,如果你想将正态分布的蒙特卡洛模拟结果进行逆转换,你可以将上述步骤反过来进行。首先,使用R矩阵和逆正交矩阵将数据反正交化;然后,使用boxcox的逆函数(使用找到的最佳lambda值)将数据反正态化;最后,使用exp函数将数据反对数化。

但请注意,这只是一个基础示例,具体的代码可能需要根据你的数据进行调整。并且,由于我无法访问到你的数据,我无法保证这个代码能够完全满足你的需求。

在Python中,我们可以这样做:

import numpy as np
import pandas as pd
from scipy.stats import boxcox
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# 生成非正态分布的时间序列矩阵
np.random.seed(42)
X = np.random.exponential(size=(100, 5))
X[:, 0] = np.random.gamma(shape=10, scale=0.1, size=100)
X[:, 1] = np.random.beta(a=2, b=5, size=100)
X[:, 3] = np.random.poisson(lam=5, size=100)
X[:, 4] = np.random.normal(loc=3, scale=2, size=100)

# 对数转换和Box-Cox转换
X_log = np.log(X + 1)
X_boxcox, _ = boxcox(X + 1)

# 正交转换(PCA)
scaler = StandardScaler()
X_standard = scaler.fit_transform(X_log)
pca = PCA()
X_pca = pca.fit_transform(X_standard)

# 蒙特卡洛模拟生成标准正态分布样本
np.random.seed(42)
mc_samples = np.random.normal(loc=0, scale=1, size=(1000, 5))

# 逆转换
mc_samples_pca = pca.inverse_transform(mc_samples)
mc_samples_log = np.exp(mc_samples_pca) - 1
mc_samples_boxcox = (mc_samples_pca * np.std(X_boxcox, axis=0)) + np.mean(X_boxcox, axis=0)

下面是将原始数据、对数转换、Box-Cox转换和正交转换后的数据进行可视化:

import matplotlib.pyplot as plt

fig, axs = plt.subplots(4, 2, figsize=(10, 10))
fig.subplots_adjust(hspace=0.5)

# 原始数据
axs[0, 0].hist(X[:, 0], bins=20)
axs[0, 0].set_title('Original data - variable 1')

axs[0, 1].hist(X[:, 1], bins=20)
axs[0, 1].set_title('Original data - variable 2')

axs[1, 0].hist(X[:, 2], bins=20)
axs[1, 0].set_title('Original data - variable 3')

axs[11].hist(X[:, 3], bins=20)
axs[1, 1].set_title(‘Original data - variable 4’)

axs[2, 0].hist(X[:, 4], bins=20)
axs[2, 0].set_title(‘Original data - variable 5’)

# 对数转换
axs[0, 0].set_ylabel(‘Log transform’)
axs[0, 0].hist(X_log[:, 0], bins=20)

axs[0, 1].hist(X_log[:, 1], bins=20)

axs[1, 0].hist(X_log[:, 2], bins=20)

axs[1, 1].hist(X_log[:, 3], bins=20)

axs[2, 0].hist(X_log[:, 4], bins=20)

# Box-Cox转换
axs[0, 0].set_ylabel(‘Box-Cox transform’)
axs[0, 0].hist(X_boxcox[:, 0], bins=20)

axs[0, 1].hist(X_boxcox[:, 1], bins=20)

axs[1, 0].hist(X_boxcox[:, 2], bins=20)

axs[1, 1].hist(X_boxcox[:, 3], bins=20)

axs[2, 0].hist(X_boxcox[:, 4], bins=20)

# 正交转换(PCA)
axs[0, 0].set_ylabel(‘PCA transform’)
axs[0, 0].hist(X_pca[:, 0], bins=20)

axs[0, 1].hist(X_pca[:, 1], bins=20)

axs[1, 0].hist(X_pca[:, 2], bins=20)

axs[1, 1].hist(X_pca[:, 3], bins=20)

axs[2, 0].hist(X_pca[:, 4], bins=20)

# 蒙特卡洛模拟结果的逆转换
axs[0, 0].set_ylabel(‘Inverse transform - Log’)
axs[0, 0].hist(mc_samples_log[:, 0], bins=20)

axs[0, 1].hist(mc_samples_log[:, 1], bins=20)

axs[1, 0].hist(mc_samples_log[:, 2], bins=20)

axs[1, 1].hist(mc_samples_log[:, 3], bins=20)

axs[2, 0].hist(mc_samples_log[:, 4], bins=20)

axs[0, 0].set_ylabel(‘Inverse transform - Box-Cox’)
axs[0, 0].hist(mc_samples_boxcox[:, 0], bins=20)

axs[0, 1].hist(mc_samples_boxcox[:, 1], bins=20)

axs[1, 0].hist(mc_samples_boxcox[:, 2], bins=20)

axs[1, 1].hist(mc_samples_boxcox[:, 3], bins=20)

axs[2, 0].hist(mc_samples_boxcox[:, 4], bins=20)

plt.show()


python实现非正态分布转正态分布
举个例子

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib
from scipy import stats
import warnings
warnings.filterwarnings("ignore")

def Box_Cox(file,sheet_name):
    df1 = pd.read_excel(file,sheet_name)
    print(df1["患者密度(人/10万人)"])
    sns.distplot(df1["患者密度(人/10万人)"],color = "#D86457")
    matplotlib.rcParams['font.sans-serif'] = ['SimHei']
    matplotlib.rcParams['axes.unicode_minus'] = False
    plt.show()

    fig = plt.figure()
    ax = fig.add_subplot(111)
    stats.boxcox_normplot(df1["患者密度(人/10万人)"], -20, 20,plot = ax)
    plt.axvline(x = stats.boxcox_normmax(df1["患者密度(人/10万人)"]),color = "#D86457")
    plt.show()

    print(stats.boxcox_normmax(df1["患者密度(人/10万人)"]))
    x = stats.boxcox(df1["患者密度(人/10万人)"],stats.boxcox_normmax(df1["患者密度(人/10万人)"]))
    sns.distplot(x,color = "#D86457")
    plt.show()

    df=pd.DataFrame(x,columns=['转换'])
    print(df)


if __name__=='__main__':
    Box_Cox("F:\医学大数据课题\论文终稿修改\实验\差异性分析.xlsx",sheet_name='人口密度分组')


Python中的numpy和scipy库,可以使用以下代码将非正态分布的归一化时间序列矩阵进行对数转换:


```python
import numpy as np
from scipy.stats import boxcox
 # 假设您的归一化时间序列矩阵为data
# 进行对数转换
log_data = np.log(data)
要进行正态转换,可以使用以下代码:


from scipy.stats import norm
 # 计算每列数据的均值和标准差
mean = np.mean(data, axis=0)
std = np.std(data, axis=0)
 # 进行正态转换
norm_data = norm.cdf((data - mean) / std)


要进行正交转换,可以使用以下代码:

from sklearn.decomposition import PCA
 # 进行主成分分析
pca = PCA()
pca.fit(data)
 # 获取主成分得分
ortho_data = pca.transform(data)

要将正态分布的蒙特卡洛模拟结果进行逆转换,可以使用以下代码:

inv_norm_data = norm.ppf(norm_data)
 # 将正交转换的得分转换回原始数据
inv_ortho_data = pca.inverse_transform(ortho_data)