从Python到R：如何在ggplot中绘制带有100条回归线的向量？

Question

我调整了一段Python代码来绘制这些图表：

这段Python代码如下：

    import pandas as pd 
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)  # For reproducibility
size = 200
R_f = 0.02  # Risk-free rate
E_R_m = 0.1  # Expected market return
beta = 1.5  # Beta for the asset
sigma = 0.05  # Standard deviation of error

# Generate market returns
E_R_m_samples = np.random.normal(E_R_m, sigma, size)

# Calculate asset returns based on the CAPM model
epsilon = np.random.normal(0, sigma, size)  # Noise
E_R_i_samples = R_f + beta * (E_R_m_samples - R_f) + epsilon  # Asset returns

# Prepare data for Bayesian regression
market_excess_return = E_R_m_samples - R_f  # Independent variable (x)
asset_excess_return = E_R_i_samples - R_f  # Dependent variable (y)

data_capm = pd.DataFrame(dict(market_excess_return=market_excess_return, asset_excess_return=asset_excess_return))

# Plot the generated data
plt.figure(figsize=(7, 7))
plt.scatter(data_capm['market_excess_return'], data_capm['asset_excess_return'], label='Sampled Data')
plt.xlabel('Market Excess Return')
plt.ylabel('Asset Excess Return')
plt.title('Generated CAPM Data')
plt.legend()
plt.show()

import pymc as pm
import arviz as az

with pm.Model() as model_capm:
    # Define priors
    sigma = pm.HalfCauchy('sigma', beta=10)
    intercept = pm.Normal('Intercept', 0, sigma=20)
    slope = pm.Normal('slope', 0, sigma=20)

    # Define likelihood
    likelihood = pm.Normal('asset_excess_return', mu=intercept + slope * market_excess_return,
                           sigma=sigma, observed=asset_excess_return)

    # Inference
    trace_capm = pm.sample(3000, return_inferencedata=True)

# Plot trace
az.plot_trace(trace_capm, figsize=(10, 7))
plt.show()

posterior_intercept = trace_capm.posterior["Intercept"].values.flatten()
posterior_slope = trace_capm.posterior["slope"].values.flatten()

# Generate a range of x values for plotting
x_plot = np.linspace(market_excess_return.min(), market_excess_return.max(), 100)

# Plot a subset of posterior predictive regression lines
plt.figure(figsize=(7, 7))
for i in range(100):  # Plot 100 posterior predictive lines
    y_plot = posterior_intercept[i] + posterior_slope[i] * x_plot
    plt.plot(x_plot, y_plot, 'c-', alpha=0.2)

# Plot observed data
plt.scatter(market_excess_return, asset_excess_return, alpha=0.5)
plt.xlabel('Market Excess Return')
plt.ylabel('Asset Excess Return')
plt.title('Posterior Predictive Regression Lines with Observed Data')
plt.show()

我用x_plot <- seq(min(market_excess_return), max(market_excess_return), length.out = 100)来模拟回归线。

不过，当我想把回归线绘制到图表中时，只能用geom_abline()添加一条拟合线，结果是这样的：

我该如何把x_plot对象包含进去呢？

我提供了完整的R代码：

library(ggplot2)
library(dplyr)
library(brms) # For Bayesian regression
library(bayesplot) # For plotting Bayesian analysis results

set.seed(42)
size <- 200
Rf <- 0.02
ERm <- 0.1
beta <- 1.5
sigma <- 0.05

ERm_samples <- rnorm(size, ERm, sigma)

epsilon <- rnorm(size, 0, sigma)
ERi_samples <- Rf + beta * (ERm_samples - Rf) + epsilon

market_excess_return <- ERm_samples - Rf
asset_excess_return <- ERi_samples - Rf

data_capm <- data.frame(market_excess_return=market_excess_return, asset_excess_return=asset_excess_return)

ggplot(data_capm, aes(x=market_excess_return, y=asset_excess_return)) +
  geom_point(aes(color='Sampled Data')) +
  labs(x='Market Excess Return', y='Asset Excess Return', title='Generated CAPM Data') +
  theme_minimal()

# Bayesian regression using brms
model_capm <- brm(asset_excess_return ~ market_excess_return, data = data_capm, 
                  family = gaussian(), prior = c(
                    prior(normal(0, 20), class = Intercept),
                    prior(normal(0, 20), class = b),
                    prior(cauchy(0, 10), class = sigma)
                  ))

# Plotting the trace
trace_capm <- as.mcmc(model_capm) 
bayesplot::mcmc_trace(trace_capm)

# Posterior predictive checks
#posterior_intercept <- fitted(model_capm)$Estimate[1,]
#posterior_slope <- fitted(model_capm)$Estimate[2,]
posterior_intercept <- fitted(model_capm)[1, "Estimate"]
posterior_slope <- fitted(model_capm)[2, "Estimate"]

x_plot <- seq(min(market_excess_return), max(market_excess_return), length.out = 100)

ggplot(data_capm, aes(x=market_excess_return, y=asset_excess_return)) +
  geom_point(alpha=0.5, color = 'tomato') +
  geom_abline(intercept=mean(posterior_intercept), 
              slope=mean(posterior_slope), 
              color="blue", 
              linetype="dashed") +
  labs(x='Market Excess Return', 
       y='Asset Excess Return', 
       title='Posterior Predictive Regression Lines with Observed Data') +
  theme_minimal()