如何用R绘制尖峰时间指数分布直方图及尖峰序列栅格图

嘿，我来帮你搞定这两个R语言的尖峰序列分析问题！

首先要明确：泊松过程生成的尖峰，其**尖峰间隔（ISI，Inter-Spike Interval）**是服从指数分布的，速率参数等于发放率fr（单位Hz）。所以我们需要先提取所有尖峰间隔，绘制直方图，再叠加理论或拟合的指数分布曲线。

步骤与代码 #

1. 从你的尖峰矩阵中提取所有试验的尖峰时间，并计算ISI：

# 提取每个试验的尖峰时间（转换为秒）
spike_times <- lapply(1:nTrials, function(trial) {
  which(SpikeMat[, trial] == 1) * dt
})

# 计算所有试验的ISI并合并
all_isi <- unlist(lapply(spike_times, function(t) {
  diff(t)  # 相邻尖峰的时间差
}))

2. 绘制ISI直方图，并叠加理论指数分布曲线：

# 绘制归一化直方图（freq=FALSE表示密度而非频数）
hist(all_isi, breaks = 30, freq = FALSE, 
     main = "尖峰间隔分布与指数拟合",
     xlab = "尖峰间隔（秒）", col = "lightblue")

# 生成理论指数分布的密度曲线
lambda <- fr  # 指数分布的速率参数等于发放率
x_seq <- seq(0, max(all_isi), length.out = 1000)
lines(x_seq, dexp(x_seq, rate = lambda), col = "red", lwd = 2)

# 添加图例
legend("topright", 
       legend = c("ISI直方图", "理论指数分布"),
       col = c("lightblue", "red"), lty = 1, lwd = 2)

如果需要用实际数据拟合指数分布（而非直接用理论值），可以借助MASS包的fitdistr函数：

library(MASS)
# 拟合指数分布
fit_result <- fitdistr(all_isi, "exponential")
# 叠加拟合的曲线
lines(x_seq, dexp(x_seq, rate = fit_result$estimate), col = "green", lwd = 2)
# 更新图例
legend("topright", 
       legend = c("ISI直方图", "理论指数分布", "拟合指数分布"),
       col = c("lightblue", "red", "green"), lty = 1, lwd = 2)

先补全你未写完的尖峰矩阵生成代码，然后用两种方法绘制栅格图：基础R和ggplot2。

第一步：补全尖峰矩阵生成代码 #

fr = 100; dt = 1/1000  # dt是毫秒转秒的系数，即0.001秒
duration = 2  # 试验时长（秒）
nBins = duration / dt  # 时间点数量，即2000
nTrials = 20  # 模拟次数

MyPoissonSpikeTrain <- function(fr=100) {
  p = runif(nBins)
  # 每个时间点发放尖峰的概率是fr*dt（发放率×时间步长）
  q = ifelse(p < fr*dt, 1, 0)
  return(q)
}

set.seed(1)
# 生成nTrials个试验的尖峰序列，转置为nBins行×nTrials列的矩阵
SpikeMat <- t(replicate(nTrials, MyPoissonSpikeTrain(fr=fr)))

方法1：基础R的image函数（快速简洁） #

# 绘制栅格图：x轴为时间，y轴为试验编号，黑色点代表尖峰
image(x = seq(0, duration, length.out = nBins), 
      y = 1:nTrials, 
      z = SpikeMat, 
      col = c("white", "black"),  # 0=白色，1=黑色
      xlab = "时间（秒）", 
      ylab = "试验编号", 
      main = "尖峰序列栅格图")

方法2：ggplot2（更美观，可自定义样式） #

需要先把宽格式的矩阵转换为长格式数据框：

library(ggplot2)
library(tidyr)

# 转换为长格式数据框
spike_df <- as.data.frame(SpikeMat) %>%
  mutate(time = seq(0, duration, length.out = nBins)) %>%
  pivot_longer(cols = -time, names_to = "trial", values_to = "spike") %>%
  # 提取试验编号（默认列名是V1、V2...，转换为数字）
  mutate(trial = as.numeric(substr(trial, 2, nchar(trial))))

# 绘制栅格图
ggplot(spike_df, aes(x = time, y = trial)) +
  geom_tile(aes(fill = factor(spike)), color = "gray") +
  # 设置颜色：无尖峰为白，有尖峰为黑
  scale_fill_manual(values = c("white", "black"), guide = "none") +
  labs(x = "时间（秒）", y = "试验编号", title = "尖峰序列栅格图") +
  theme_minimal()

内容的提问来源于stack exchange，提问作者14thTimeLord