This study evaluates paleoclimate sensitivity over the past 800,000 years from proxy-based reconstructions of changes in global temperature, ice sheets and sea level, vegetation, dust, and greenhouse gases. This analysis uses statistical methods that are not biased by the variable (heteroscedastic) uncertainty in the reconstructions, and applies a Monte Carlo-style probabilistic framework to quantify several sources of measurement and structural uncertainty. Not addressing the heteroscedastic uncertainty would result in regression results that underestimate paleoclimate sensitivity by over 30%, and not using a probabilistic framework could underestimate the credible interval by fivefold. A comparison of changes in global temperature (ΔT) and changes in radiative forcing from greenhouse gases, ice sheets, dust, and vegetation (ΔR[GHG,LI,AE,VG]) over the past 800 kyr finds that the two are closely coupled across glacial cycles with a correlation of 0.81 (0.6 to 0.9, 95% credible interval). The variation of ΔT with ΔR over the past 800 kyr is non-linear, with lower correlation and lower responsiveness at colder temperatures. The paleoclimate sensitivity parameter estimates (S[GHG,LI,AE,VG]) are 0.84 °C/W/m2 (0.20 to 1.9 °C/W/m2, 95% interval) for interglacial periods and intermediate glacial climates and 0.53 °C/W/m2 (0.08 to 1.5 °C/W/m2, 95% interval) for full glacial climates, 37% lower at the median. The estimates of S[GHG,LI,AE,VG] and the pattern of state dependence are similar across glacial cycles over the past 800 kyr. This analysis explicitly includes several sources of uncertainty and is still able to provide a strong upper bound for the paleoclimate sensitivity parameter for interglacial periods and intermediate glacial climates: over 1.5 °C/W/m2 is < 10% probability, 1.7 °C/W/m2 is < 5% probability, and over 1.9 °C/W/m2 is < 2.5% probability.